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Abstract 

This report was written in 2001 and it is a translation of work that was originally published in Greek, in the 
author's diploma thesis in July 1998. It may contain minor mistakes and should not be considered a complete 
study. It however touches upon several of the considerations that will be included in the complete paper. The 
work was done jointly with E. Gallopoulos and was presented at the FOCM99 conference, at the 5th IMACS 
conference on iterative methods in scientific computing, and at the 50th annual meeting of SIAM. 

We present an iteration for the computation of simple eigenvalues using a pseudospectrum approach. The most 
appealing characteristic of the proposed iteration is that it reduces the computation of a single eigenvalue to a small 
number of eigenvalue computations on Hermitian matrices. We show that this number is directly associated with the 
matrix pseudospectrum. We present numerical results and we discuss advantages and drawbacks of the method. We 
also discuss its relationship with an iteration that was proposed independently in [Stewart, O' Leary, ETNA, Vol 8, 
1998]. 

1 A pseudospectrum setting 

We shortly describe some well known characteristics of matrix pseudospectrum. Denote with A(A) the spectrum of a 
matrix A, with A e (A) the e-pseudospectrum of A, with cr m in = &n < cjv-i • ■ • < f l the singular values of A, and 
with D(z, q) and D°(z, g) a closed disk and an open disk respectively, with center z and radius g. Let N denote a 
normal matrix. 

Theorem 1 

• A(A) C A t {A)fore > 0. 

• e < ei 4^ A £ (A) C A ei (A). 

• A(N) = A(A) A e (N) C A e (A). 

These properties state that e-pseudospectrum forms closed curves around eigenvalues. The crucial property of 
subharmonicity of the norm of the resolvent \ \(zl — A)^ 1 1| (see J2J), assures that 

Vmin(zl — A) — local minimum = 0«z£ A(A) (1) 

This alternative characterization of an eigenvalue together with the properties of matrix pseudospectrum appeal for the 
corresponding optimization problem: 

min/(z) = a min (zl -A), z E C 

Notice that in this formulation we are seeking for a minimizer over the complex plane. This is not usually the 
case in other methods for the computation of eigenvalues . The "heart" of any competent minimization algorithm can 
naturally be the following theorem by Sun ifTTl . 
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Theorem 2 Denote with <Ji the i singular value and Ui, Vi the corresponding left and right singular vectors. Let 
also z = x + iy G C\ A(A). Then g(x, y) — a^zl — A) is real analytic in a neighborhood of '{x , y), if a m in(zl — A) 
is a simple singular value. For the gradient of g(x, y) we have 

Vs(x,y) = (Re(vtui),Im(v*Ui)) 

What is interesting, is that Vcr m i n (2;/ — A) can be computed almost as easy as er m i n itself. Given that cr m i n can be 
computed by means of some of the well known methods for Hermitian matrices, we can readily setup the framework 
for a minimization algorithm. Using for example a steepest descent strategy we get the extremely compact Mat lab 
routine shown below 



Eigenvalue via Hermitian Computations 

function z=sdeig (A, starting_point) 

z= starting_point ; s_min= 1; 
while (s_min > tol) 

[u_min, s_min, v_min] = svds (z*I-A) ; 

direction = vjnin *u_min / |v_min *u_min|; 

step_size = ? ; 

z = z - st ep_si ze *direct ion ; 

end 



The determination of the step size is obviously of critical importance for the effectiveness of the sdeig. In our 
experiments we used the general purpose minimization function fminu of Matlab Optimization Toolbox, 
supplemented with the gradient information. Alas, even for a moderate accuracy of the output approximate eigenvalue, 
the number of iterations (and (J mnl evaluations) of fminu, usually exceeded the number of 100. 

The main part of this paper aims to determine a strategy for selecting the step size. Using mainly results from our 
work for the efficient computation of the pseudospectrum [3], we give a simple and computationally efficient formula 
for the step size, which as we shall see, when incorporated in sdeig gives a powerful iteration. 

2 Determining Step Size 

In this section we will use theorems proven in [3 1. Denote by D° (z, r) the open disk with center z and radious r. We 
usually denote with a m i n (z) the minimum singular value of zl — A and with u m i n (z), v ni i n (z) the corresponding left 
and right singular vectors. 

Theorem 3 If a min (zl — A) = r > e, then D°(z,r — e) R A e (A) = 0. In the extreme case e — 0, we have 
D°(z, r) n A(A) = 0. 

Is fj m in (zl — A) a good candidate for step size? The answer is affirmative in the case of normal matrices, as the 
following corollary shows: 

Corrolary 1 If N 6 C I1X " is normal, Zi,i — 1,. .. ,n are its eigenvalues, and o~ m i n (zI — N) = r > e then the 
perimeter dD(z, r — e) contains at least one point of A e (N). In the extreme case e = 0, the disk D(z,r) "touches " at 
least one eigenvalue of N. 

Corollary Q] along with the observation that pseudospectra of normal matrices are concentric circles around eigen- 
values, guarantee that sdeig always finds in one step the eigenvalue closest to the starting point. However, in 
non-normal matrices the pseudospectra expansion can be very fast. In other words, the value of er m i n (z) can be very 
small comparing to the distance of z from A(A), so that the convergence of sdeig is very slow. Our basic effort will 
be to measure and cancel this fast pseudospectra expansion. To this aim we define the pseudospectrum sensitivity at 
a point z to be 
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Definition 1 

pss(z) = \v" m (z) * UminO)!" 1 (2) 

Using the above definition the effectiveness of Theorem[3]can be extended as follows. 

Theorem 4 Let z — i + 6 C \ A(A) and L the line segment between z and the closest eigenvalue. Ifpss(z) < 
pss(z') for every point z' £ L, then D°(z, a m - ul (z) ■ pss(z)) n A e (A) = 0. 

Proof Let g(x, y) = cr min ((x + yi)I - A). By definition g(x, y) - g(x ,y ) = ^ Jj* o) Vg(x, y)n L dL , where n L 
is a normalized vector parallel to L. Then using Theorem|2] we have 

r{xi,Vi) r(xi,vi) 

\g(xo,yo) -g(xi,yi)\ = \ & Vg(x,y) ■ n L dL\ < j> \V 'g(x,y) ■ n L \ dL 
pss(z)^ 1 dL < pss(zo)^ 1 (p dL = pss(zo)^ 1 ■ \zq — Zi\. 



Of course, the key question concerning [4] is its applicability, given that its assumption does not hold in general. 
Notice first that even in cases where the assumption does not hold, the statement of the theorem can be true, since 
pss(z) is not in general constant in L. In (3) we showed that the above statement is in general true, with the exception 
of some areas between eigenvalues. We could briefly explain this phenomenon by means of some basic theorems and 
observations. 

Theorem 5 Let zq be a simple eigenvalue of A with y, x the corresponding left and right eigenvectors. If we define 
the eigenvalue condition as k{zq) — ^T-4^jrp^ , then 

lim pssa{z) = k(zq)^ 1 

z^zo 

Observe now that in the limit \z\ — > +oo, we have v^u n — > 1. This fact along with Theorem[5]and continuity of 
pss(z) is an indication of the pss(z) behavior with dist(z, A(A)) increasing. Of course, our claim gives no guarantee 
for the pss(z) behavior close to A(A). The following theorem taken from [10] gives a local description around simple 
eigenvalues. 

Theorem 6 Let z e be a simple eigenvalue of A with y, x the corresponding left and right eigenvectors and A — A + E 
be a perturbation of A. Then there is a unique eigenvalue z e of A such that 

\v A 

Second order terms in Theorem [6] translate directly in a decrease of pss(z) locally around simple eigenvalues. 
Concerning areas between eigenvalues where the statement of Theorem|4]fails, we can use a basic result from [[8). 

w min( z ) * u xmn{z) = 3E such that z is a double eigenvalue of A + E. (4) 

In that case a n — a n -\ is a necessary condition. Since pss(z) is continuous there are areas around z where pss(z) 
is very small and thus gives very large incorrect exclusions. We will defer further discussion on the behavior of 
pseudospectrum sensitivity behavior until section [7] The final argument of this section is given in Figure Q] which 
essentially shows that the assumption of Theorem|4]is very weak. 
We have set grounds which justify our selection for the step size: 

step.size — a mlu (z) ■ pss(z) (5) 
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grcar(32} pseudo spectrum sensitivity 



pentoep(3Z,0.Q.5,0,0,1) pseudDspccTum iivity 




Figure 1: Contour plots of pss(z) for matrices grcar(32) and propeller(32) . 

3 Numerical Results 

In this section we present numerical results that in principal show the effectiveness of sdeig. Starting points have 
been selected using the enclose algorithm of [3|. We use the complete svd decomposition, but essentially the 
same results are obtained with Arnoldi methods, for example with the function svds of Mat lab. In the following, 
with A we denote an eigenvalue computed by eig, with "eig_sens" we denote the sensitivity of the corresponding 
eigenvalue computed by eigsens. All test matrices and function eigsens are from the Testmatrix Toolbox 
for Mat lab. 

The steepest descent is stopped when the first increase of the numerical values of cr m j n occurs. In that sense, 
the implementation of a stopping criterion is straightforward. However as we can see in Table Q] convergence may 
have been obtained earlier. Hence, a better implementation of a stopping criterion could stop avoid some of the last 
iterations. An example of the steepest descent process and a zoom in the area of the eigenvalue is given in Figure|2] 



Trace of the Steepest Descent, Matrix pentoep(32.0.( 



te Steepesl Descent, Matrix pentoep(32,0, 



Figure 2: The steepest descent process 



In Table|2], we see that the number of iterations can be decreased with the choice of a "good" starting point (15/8), 
and of an appropriate stopping criterion tol (17/12). 

In the examples of Table[3], we use a slightly different matrix changing a single parameter in function pentoep. 
The spectrum of the new matrix is of similar form with that discussed on Table [2] but eigenvalue sensitivities are 
clearly reduced. We use the starting points of Table|2] The effects are clear in the size of |A — A| and more importantly 
in the number of iterations. 
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Starting Point 


Iterations 


A 


|A-A| 


s_min 


eigjsens 


0.6+0.5i 


15 


6.96686e-01 


3.84137e-014 


4.933e-018 


9.45e+003 


0.6+0.5i 


14 




9.23759e-014 


1.383e-017 




0.6+0.5i 


12 




1.29933e-006 


1.374e-010 





Table 1: Matrix pentoep(32,0,0.5,0,0,l) 



Starting Point 


Iterations 


A 


|A-A| 


s_min 


eigjsens 


1.0+O.li 


8 


1.02409e+00 


1.90412e-014 


8.191e-019 


2.75e+004 


1.0+0.5i 


15 


9.32779e-01 


1.85407e-014 


2.091e-019 


2.32e+004 


1.3+0.0i 


17 


1.17998e+00 


5.10702e-015 


6.897e-019 


5.21e+003 


1.3+0.0i 


12 




1.35447e-014 


1.770e-018 




-0.4+0.5i 


9 


-(3.48-6.03i)e-01 


3.57989e-014 


2.018e-018 


9.45e+003 



Table 2: Matrix pentoep(32,0,0.5,0,0,l) 



Starting Point 


Iterations 


A 


|A-A| 


s_min 


eigjsens 


1.0+O.li 


7 


1.03032e+00 


4.04639e-022 


1.586e-018 


2.82e+002 


1.0+0.5i 


12 




2.22044e-016 


4.832e-019 




-0.4+0.5i 


7 


-(3.44-6.00i)e-01 


3.51803e-016 


3.08e-018 


1.86e+002 



Table 3: Matrix pentoep(32,0,0.7,0,0,l) 



Starting Point 


Iterations 


A 


|A-A| 


s_min 


eigjsens 


1.0+0.2i 


10 


9.32039e-01 


8.88178e-015 


2.727e-017 


9.09e+001 


0.8+0.2i 


13 


7.03348e-01 


1.14330e-012 


3.041e-017 


3.98e+003 


0.6+0.2i 


15 


4.94699e-01 


8.37324e-012 


2.023e-017 


7.38e+004 


0.4+0.2i 


17 


3.47946e-01 


3.01335e-011 


1.727e-017 


1.08e+006 


0.2+0.2i 


18 


1.72127e-01 


1.71558e-011 


1.634e-019 


3.26e+007 



Table 4: Matrix kahan(32) 
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In the examples of Table [4] we use a nice property of matrix kahan(32). Its eigenvalues lie on the real axis in the 
range [0.11, 1.0.1], and their sensitivities satisfy ei < ej <^> sens(ei) > sens(ej). We observe that the increase of 
eigenvalue (and pseudospectrum) sensitivity reflects in a monotone way in the number of iterations and in the size of 
| A — A | . However, an O(10 6 ) increase of the eigenvalue sensitivity does not even double the number of iterations. We 
also observe that independently from eigenvalue sensitivity, a m i n (z — A) < 10~ 16 . For a variety of similar results 
refer to [?]. 



4 The role of clustering 

In the previous section we observed a connection between the number of iterations and eigenvalue sensitivity (or 
closeness of the original matrix to defective matrices - see an argument in iflOl ). But, how clusters of eigenvalues on 
the complex plane can affect the number of iterations? We settle the question by means of a numerical example, and 
an extreme case theorem. 



Spectrum ol matrices rattd{t00| and rand(32J 
3| 1 1 1 1 r- 




Figure 3: Examples of geometric clusterings 

In Figure[3]we plot spectra of rand ( 32 ) and rand ( 1 ) . It is obvious that eigenvalues of rand ( 100 ) exhibit 
a more intense geometric clustering. However, sensitivity of eigenvalues is less than 100 in both cases. In particular 
mean (eigsens (rand (32) ) ) =4 . 59 and mean (eigsens (rand (100) ) ) =6 . 25. For a small number of 
starting points inside the cluster, we obtained mean ( iterations ( rand ( 32 ) ) ) =6. 8 and 
mean (iterations ( rand (100) ) )=6.9. We now give an extreme case theorem. 

Theorem 7 Let N be a normal ' matrix with eigenvalues Ai, .., X n , with Ai| > .. > A„|, and singular values a\ > 
.. > a n . The following hold: 

1. | A, ; | = a. t . 

2. For every z £ C, if zl — N = UYjV h , there is a permutation tt, such that 

=z- S(i,i) * V(;,i)' *U(: A), i = 1, ...,n. 

Proof For claim 1, see 0. In the case of normal matrices we have u^Vi for every singular value a^. Let Zi be 
the eigenvalue for which \zi — z\ = Oi{z). Let L be the line between z and Z{. For each point z' of L, there is some 
j < i such that \z' — Zi\ = <Jj(z!l — A). We now split L in segments where the index j is constant and we consider 
the integral, as we did in proof of Theorem |4] for the corresponding <jj along each segment. Observe that since N is 
normal, we have Ui(z') H Vi(z') = 1, for every i, z' ', and the inequalities of the proof become strict equalities. Hence, 
the gradient of <jj must be parallel to L at every point of it. □ 
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computed eigenvalues 


\^eig ~ M 


4.354989548789193e+000+2.213138973009386e-001i 
3 .9024 1 77 14058627e+000 +3.6 1 9844954770 19 1 e+OOOi 
2.750755557087913e+000+9.117844570868215e+000i 
3 .34486 1 206758000e+000 +9.20404 1 639930084e+000i 
7.1368 1 594476 1 636e+000 +8 .562489580579660e+000i 


1.766e-015 
3.203e-015 
1.432e-014 
7.105e-015 
1.387e-014 



Table 5: Eigenvalues by means of svd 



In Table [4] we give the eigenvalues of a normal 5x5 matrix with no trivial eigenvectors, computed using the 
algorithm of Theorem|7] for the point z = 1.5. We also compare with values computed by eig. 

A well known result (see [4|), is that svd requires (by a factor) less flops than eig. In our 5x5 example, 
flops (eig) = 4695, flops (svd) =2293. Theorem [7] directly suggests a way for exploiting normality in the 
computation of eigenvalues of normal matrices with direct or indirect methods. 

This section established also, that the size of the matrix does not affect the number of iterations. 

5 Using Exclusions 

In the previous sections we established sdeig, a steepest descent algorithm for the computation of a simple eigen- 
value. The central observation of this section is that the steepest descent procedure can give additional useful informa- 
tion. We repeat the key claim of Section|2] 

Claim 1 D°(z, cr m i n (zl — A) ■ pss(z)) and A(A) have under weak conditions no common points. 

This claim implies that we can use all the cr m - ln {z) computations of sdeig to obtain regions of the complex plan 
which do not contain eigenvalues. Figure [4] shows the area excluded after two runs of eig. The use of "localization 
information" could be an additional feature of an eigenvalue algorithm which will use our ideas. Of course, exclusions 
will be most informative in the case of hermitian and in general, normal matrices. 



neig and exiusions 




-1 -0.5 0.5 1 1.5 

real 



Figure 4: Using Exclusions 

In Figure [4] also observe that exclusion disks intersect. A more sophisticated algorithm could choose a slightly 
bigger step size (at least in areas away from eigenvalues) in order to diminish the number of iterations. A natural 
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way to try such an improvement would be to use second order information for <7 m i n . These are open issues for future 
research. 



6 The Oleary-Stewart Result 

The iteration of sdeig is given in |9| slightly disguised. It is interesting that the iteration is derived from a very 
different perspective. The starting point is the following variant of the Rayleigh quotient method. Let A be a simple 
eigenvalue of A with right and left eigenvectors x and y H . Let v and w H be approximations to x and y H , and let r be 
an approximation to A. The new approximations v, w H and f are generated as follows: 

1. v = {A-riyH 

2. w H = w H (A - tI)- 1 

3. f = w H Av/w H v 

The precedure is then iterated. The basic observation of (9J is that if r is close to an eigenvalue the singular vectors 
should approximate the eigenvectors. Using only this observation the Rayleigh quotient iteration is transformed to the 
iteration of sdeig. Finally, it is given a proof that the iteration achieves local quadratic convergence. 

Theorem 8 ife = |A — r| is sufficiently small, there is a constant fi such that: 

e = |A - f| < C(( £ / M ) 2 + (e/ M )(||<MI + ||<U) + , (6) 

where 

- , ^ (7, 

VI - (t/M) 2 - ls(2(e//i) + + Ifcll) 

The constant /i is the lower bound of a n -\{z) in a neighborhood of \ and 5 U , S w denote inaccuracies in the computa- 
tion of the singular vectors u, w. 

Unfortunately, the result of [9| holds only for (unbounded) small neighborhoods around simple eigenvalues, and 
fails to explain the general convergence exhibited by sdeig. It is however useful in proving rigorously the ultimate 
convergence of sdeig, and establishing a better stopping criterion. 



7 Pseudospectrum Sensitivity Behavior 

Figure Q] shows that pseudospectrum sensitivity exhibits a behavior more structured than that implied from theorems 
and observations in Section|2] In this section we will attempt to give a more accurate picture of that structure. We first 
give a general theorem taken from [1] (see also ||6l). 

Theorem 9 Let X be a Banach space and Q an open subset ofC. Consider a function f : f2 — > X holomorphic in f2 
and continuous in Q. Then \ \f\ \ is subharmonic in fl 

We will need one more basic result contained in |7| (see page 84). 

Theorem 10 If A(z) £ c«x™ { s an holomorphic function of z £ C, and X(z) is a simple eigenvalue of A(z) in some 
domain C C then X(z) is holomorphic in f2. Moreover, the corresponding left and right eigenvectors y H (z), x are 
holomorphic in fl. 

We are now ready to give the main theorem. 

Theorem 11 If A G C" x ™, pss(z) is subharmonic in every domain £1 £ C, where cr m i n (z) is simple. 
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Proof Let zl— A = U(z) H Y,{z)V{z). It is well known that the hermitian matrix (zI—A) H (zI—a) has eigenvalues 
&i(z) > ■ ■ ■ > c„(z), and corresponding left and right eigenvectors u±(z), . . . , u n (z) and Vi(z), . . . , v n (z). Clearly 
a n (z) = cr m i n (z) is simple if u n -i(z) ^ cr n (z). In that case, from TheoremfTUl <j n and Ui(z) H ,Vi(z) are holomorphic 
in z . Hence uf(z) ■ Vi(z) is holomorphic. Subharmonicity follows from Theorem[9]and the definition of pss(z). □ 

The Maximum Principle is an elementary property of subharmonic functions. Specifically, the following theorem 
holds. 

Theorem 12 Let be a bounded domain in C and dfl be its boundary. Suppose f is subharmonic in Q. Then 

sup/(s) = max/(z) (8) 
sen CeSfi 

The Maximum Principle can be directly applied to pss(z). Notice also, that since pss(z) = =>• <r„ = cr„_i, 
the function pss(z)^ 1 is subharmonic in the same domain. Hence, the Maximum Principle can also be applied to 
pss(z)~ 1 . These observations are summarized in the following corollary. 

Corrolary 2 Let Q = C \ {z : a n (z) ^ a n -i(z)}. Then 

sup pss(s) = max pss(z) 
sen cedn 

inf pss(s) = min pss(z) (9) 
sen cean 

We can strengthen the result by taking into consideration the behavior of pss(z) close to points where a n (z) = 
tj n -i(z). These are exactly the points where disjoint curves of A e (A) collide. 

Corollary |2] together with some of the observations of Section |2] plausibly explains why pss(z) decreases in a 
monotonic way as dist(z, A(A)) increases. However, the similarity of pss(z) contours with A e contours is striking 
and obviously more structured. We will need the extension of Sun's Theorem. 

Theorem 13 Let U(x,y)T,V H (x,y) = A — (x + iy)I, where a — cr n (xo,yo) is simple and non-zero. Let U = 
(u\U2 ■ ■ ■ V = (t>ii>2 ■ ■ ■ « n _i) and E = diag(cri, . . . , er n _i). We then have 

^(x Q ,y Q ) = Re{r H <S>r + l H <5>l + 2l H ^r}+Im{u%v n } 2 /a ni 
d 2 a 

dxdy ( Xo,Vo ^ = - Im {2' lH ^ r } + Im i u n v n }Re{u" v n }/a n , 
d 2 a 

— (x ,yo) = Re{r H ^r + l H ^l-2l H ^r} + Re{u^v n } 2 /a n , 

where $ = (T„(cr^J - E 2 ) -1 ,* = t{a 2 J - E 2 )" 1 , r = U H v n and I = V H u n . 

Unlike the first order case, computation of second order derivatives is prohibitive for practical purposes. However, 
it gives a very good tool for further theoretical investigation. In Figure [5] we plot the determinant of the Hessian of 

CminM. 

Obviously, the pattern is similar with that of A e and pss(z), but clearly, we cannot use a subharmonicity argument. 
However, we can observe that expansion is slower than that of pss(z). Intuitively, we could expect that the sensitivity 
of analogous constructions for derivatives of order k asymptotically approaches 1 when k — > oo. A quantity that can 
naturally capture details about these expansion phenomena is the pseudospectrum area, defined in the obvious way. 

In Figure|6]we measure the area (as the number of grid points) of A c as a function of e, of matrix pentoep(32, 0, a, 0, 0, 1), 
for three different values of parameter a. This family of matrices has the nice property that the spectra are almost 
identical for different values of a, but the sensitivities of eigenvalues change significantly. Thus, the A c area depends 
mainly on pss(z). We also measure the area of A e for the matrix grcar(32). It is clear that a general "inverse expo- 
nential law" shows up, in the sense that despite the initial explosion, the expansion of pseudospectra stabilizes quickly. 
The curves are slightly more steep for bigger sensitivities of eigenvalues. 

We now give some theorems which give a way for reasoning about the pseudospectra expansion and its relation to 
pss(z). First, a general theorem of Konrod. 
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Figure 5: Contour plots of det(H(<r m i n (z))) for matrices grcar(32) and propeller(32) . 



pentoep(32,0,o,0,0,1), 




grcar(32). pseud ospectrun 



0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 



0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Figure 6: Evolution of pseudospectrum area 



Theorem 14 Let G C IR. n be a domain and f : G — > R a smooth function, E t — {x G G\f(x) = t} the level set of 
the function f and ds the in — 1)- dimensional surface element on E t . Then 



ds 



/■max / 

meas{G) =L f { j Et iv/i 



)dt 



Theorem[T4]directly applies to our case giving the following corollary. 
Corrolary 3 Let dA t be the boundary of A t and dz be the 1 -dimensional surface element on dA t . Then 

area(A e (A)) = / ( / {pss(z)dz))dt 

JO JdA t 
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